This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

K-mean

head(iris)
plot(iris[,3:4])


irisCluster <- kmeans(iris[, 3:4], 3, nstart = 20)
irisCluster$cluster <- as.factor(irisCluster$cluster)
library(ggplot2)

ggplot(iris, aes(Petal.Length, Petal.Width, col = irisCluster$cluster, fill=Species)) + 
  stat_ellipse(geom = "polygon", col = "black", alpha = 0.5) +
  geom_point()


k.max <- 15
wss <- sapply(1:k.max, function(k){kmeans(iris[, 3:4], k, nstart=50, iter.max = 15)$tot.withinss})

plot(1:k.max, wss,
     type="b", pch = 19, frame = FALSE, 
     xlab="Number of clusters K",
     ylab="Total within-clusters sum of squares")

Periodogram

counter <- 500
x <- seq(counter)

y <- 2 * cos(2 * pi * (1/50) * x + 0.6 * pi)
signal <- y + rnorm(counter)
plot(x, y, type = "l")


plot(x, signal, type = "l")

sunspots <- scan("/Users/ansonleung/Desktop/Applied Financial Engineering/Financial Times Series Analysis/Lecture 3/sunspots.dat")
Read 459 items
plot(sunspots, type="b")

x = diff(sunspots)

I = abs(fft(x)/sqrt(458))^2
P = (4/458)*I[1:230]

freq = (0:229)/458

plot(freq, P, type="l")


library(astsa)

x <- scan("/Users/ansonleung/Desktop/Applied Financial Engineering/Financial Times Series Analysis/Lecture 3/recruit.dat")
Read 453 items
mvspec(x, log="no")


k = kernel("daniell", 4)
mvspec(x, k, log="no")


k = kernel("daniell", c(4,4)) 
mvspec(x, k, log="no")



specvalues = mvspec(x, k, log="no")

specvalues$details
       frequency   period  spectrum
  [1,]    0.0021 480.0000 3383.5058
  [2,]    0.0042 240.0000 3723.6683
  [3,]    0.0063 160.0000 4063.3578
  [4,]    0.0083 120.0000 4498.0099
  [5,]    0.0104  96.0000 4948.5148
  [6,]    0.0125  80.0000 5497.4841
  [7,]    0.0146  68.5714 5948.6416
  [8,]    0.0167  60.0000 6450.8076
  [9,]    0.0187  53.3333 6840.4396
 [10,]    0.0208  48.0000 7154.9638
 [11,]    0.0229  43.6364 7198.2791
 [12,]    0.0250  40.0000 7093.4142
 [13,]    0.0271  36.9231 6889.4382
 [14,]    0.0292  34.2857 6761.5608
 [15,]    0.0312  32.0000 6389.7518
 [16,]    0.0333  30.0000 6060.4244
 [17,]    0.0354  28.2353 5720.6943
 [18,]    0.0375  26.6667 5252.6571
 [19,]    0.0396  25.2632 4727.8013
 [20,]    0.0417  24.0000 4252.4430
 [21,]    0.0437  22.8571 3848.0378
 [22,]    0.0458  21.8182 3493.5831
 [23,]    0.0479  20.8696 2963.9326
 [24,]    0.0500  20.0000 2572.4034
 [25,]    0.0521  19.2000 2161.3083
 [26,]    0.0542  18.4615 1700.2617
 [27,]    0.0562  17.7778 1389.7388
 [28,]    0.0583  17.1429 1146.9481
 [29,]    0.0604  16.5517  964.3456
 [30,]    0.0625  16.0000  823.3735
 [31,]    0.0646  15.4839  686.7090
 [32,]    0.0667  15.0000  897.8008
 [33,]    0.0688  14.5455 1103.4691
 [34,]    0.0708  14.1176 1328.1553
 [35,]    0.0729  13.7143 1607.2885
 [36,]    0.0750  13.3333 1837.4561
 [37,]    0.0771  12.9730 2074.8114
 [38,]    0.0792  12.6316 2302.5845
 [39,]    0.0813  12.3077 2516.8509
 [40,]    0.0833  12.0000 2746.4671
 [41,]    0.0854  11.7073 2461.0304
 [42,]    0.0875  11.4286 2177.9664
 [43,]    0.0896  11.1628 1898.9515
 [44,]    0.0917  10.9091 1623.8459
 [45,]    0.0938  10.6667 1384.6294
 [46,]    0.0958  10.4348 1137.8088
 [47,]    0.0979  10.2128  899.7337
 [48,]    0.1000  10.0000  661.3318
 [49,]    0.1021   9.7959  398.7428
 [50,]    0.1042   9.6000  401.7536
 [51,]    0.1063   9.4118  402.7817
 [52,]    0.1083   9.2308  401.9023
 [53,]    0.1104   9.0566  392.1142
 [54,]    0.1125   8.8889  359.3417
 [55,]    0.1146   8.7273  324.0378
 [56,]    0.1167   8.5714  289.3280
 [57,]    0.1188   8.4211  258.3219
 [58,]    0.1208   8.2759  242.3243
 [59,]    0.1229   8.1356  226.1004
 [60,]    0.1250   8.0000  214.9562
 [61,]    0.1271   7.8689  210.3355
 [62,]    0.1292   7.7419  211.2108
 [63,]    0.1312   7.6190  224.6399
 [64,]    0.1333   7.5000  237.9108
 [65,]    0.1354   7.3846  250.1213
 [66,]    0.1375   7.2727  265.0379
 [67,]    0.1396   7.1642  270.3521
 [68,]    0.1417   7.0588  261.7040
 [69,]    0.1437   6.9565  244.2570
 [70,]    0.1458   6.8571  218.1679
 [71,]    0.1479   6.7606  190.8703
 [72,]    0.1500   6.6667  183.1164
 [73,]    0.1521   6.5753  180.4624
 [74,]    0.1542   6.4865  182.3508
 [75,]    0.1562   6.4000  186.0763
 [76,]    0.1583   6.3158  195.4318
 [77,]    0.1604   6.2338  214.3466
 [78,]    0.1625   6.1538  240.0628
 [79,]    0.1646   6.0759  267.9256
 [80,]    0.1667   6.0000  296.1027
 [81,]    0.1687   5.9259  284.9247
 [82,]    0.1708   5.8537  270.5449
 [83,]    0.1729   5.7831  248.7994
 [84,]    0.1750   5.7143  221.6484
 [85,]    0.1771   5.6471  195.8172
 [86,]    0.1792   5.5814  170.1639
 [87,]    0.1812   5.5172  141.6886
 [88,]    0.1833   5.4545  114.6808
 [89,]    0.1854   5.3933   86.7627
 [90,]    0.1875   5.3333   79.8203
 [91,]    0.1896   5.2747   73.7594
 [92,]    0.1917   5.2174   70.9909
 [93,]    0.1937   5.1613   68.8111
 [94,]    0.1958   5.1064   63.6917
 [95,]    0.1979   5.0526   58.3660
 [96,]    0.2000   5.0000   54.4895
 [97,]    0.2021   4.9485   50.4913
 [98,]    0.2042   4.8980   47.7787
 [99,]    0.2062   4.8485   43.3229
[100,]    0.2083   4.8000   39.8791
[101,]    0.2104   4.7525   36.2312
[102,]    0.2125   4.7059   34.4717
[103,]    0.2146   4.6602   34.1750
[104,]    0.2167   4.6154   35.9153
[105,]    0.2188   4.5714   37.6983
[106,]    0.2208   4.5283   39.0949
[107,]    0.2229   4.4860   40.4487
[108,]    0.2250   4.4444   41.2742
[109,]    0.2271   4.4037   41.9460
[110,]    0.2292   4.3636   43.1031
[111,]    0.2312   4.3243   44.1189
[112,]    0.2333   4.2857   48.0360
[113,]    0.2354   4.2478   48.5140
[114,]    0.2375   4.2105   49.8358
[115,]    0.2396   4.1739   51.4480
[116,]    0.2417   4.1379   52.1465
[117,]    0.2438   4.1026   53.8571
[118,]    0.2458   4.0678   54.2340
[119,]    0.2479   4.0336   54.5173
[120,]    0.2500   4.0000   54.5018
[121,]    0.2521   3.9669   49.3128
[122,]    0.2542   3.9344   46.0916
[123,]    0.2562   3.9024   41.1002
[124,]    0.2583   3.8710   36.2660
[125,]    0.2604   3.8400   31.9392
[126,]    0.2625   3.8095   27.6790
[127,]    0.2646   3.7795   24.6577
[128,]    0.2667   3.7500   21.4706
[129,]    0.2688   3.7209   18.9509
[130,]    0.2708   3.6923   18.8600
[131,]    0.2729   3.6641   18.8698
[132,]    0.2750   3.6364   19.9140
[133,]    0.2771   3.6090   21.0819
[134,]    0.2792   3.5821   23.3498
[135,]    0.2812   3.5556   25.4372
[136,]    0.2833   3.5294   27.7179
[137,]    0.2854   3.5036   30.2683
[138,]    0.2875   3.4783   32.2368
[139,]    0.2896   3.4532   33.8993
[140,]    0.2917   3.4286   34.4704
[141,]    0.2938   3.4043   34.6390
[142,]    0.2958   3.3803   34.5751
[143,]    0.2979   3.3566   32.5440
[144,]    0.3000   3.3333   30.7099
[145,]    0.3021   3.3103   28.2449
[146,]    0.3042   3.2877   26.0939
[147,]    0.3062   3.2653   24.6177
[148,]    0.3083   3.2432   23.7187
[149,]    0.3104   3.2215   23.9083
[150,]    0.3125   3.2000   24.7318
[151,]    0.3146   3.1788   25.6631
[152,]    0.3167   3.1579   27.6031
[153,]    0.3188   3.1373   29.7384
[154,]    0.3208   3.1169   31.8160
[155,]    0.3229   3.0968   32.9871
[156,]    0.3250   3.0769   34.3741
[157,]    0.3271   3.0573   35.3200
[158,]    0.3292   3.0380   35.5154
[159,]    0.3312   3.0189   34.8757
[160,]    0.3333   3.0000   34.1192
[161,]    0.3354   2.9814   33.1808
[162,]    0.3375   2.9630   31.5040
[163,]    0.3396   2.9448   29.6633
[164,]    0.3417   2.9268   28.9341
[165,]    0.3438   2.9091   27.0342
[166,]    0.3458   2.8916   25.3162
[167,]    0.3479   2.8743   23.8679
[168,]    0.3500   2.8571   22.9800
[169,]    0.3521   2.8402   22.1128
[170,]    0.3542   2.8235   20.9879
[171,]    0.3562   2.8070   20.8163
[172,]    0.3583   2.7907   21.1241
[173,]    0.3604   2.7746   20.4429
[174,]    0.3625   2.7586   20.1445
[175,]    0.3646   2.7429   20.2149
[176,]    0.3667   2.7273   20.4241
[177,]    0.3687   2.7119   20.7374
[178,]    0.3708   2.6966   21.0242
[179,]    0.3729   2.6816   21.4956
[180,]    0.3750   2.6667   21.2115
[181,]    0.3771   2.6519   20.8832
[182,]    0.3792   2.6374   20.5005
[183,]    0.3812   2.6230   20.2490
[184,]    0.3833   2.6087   19.6282
[185,]    0.3854   2.5946   19.0138
[186,]    0.3875   2.5806   18.2154
[187,]    0.3896   2.5668   17.5413
[188,]    0.3917   2.5532   17.2190
[189,]    0.3938   2.5397   17.0795
[190,]    0.3958   2.5263   16.8594
[191,]    0.3979   2.5131   17.0401
[192,]    0.4000   2.5000   17.3831
[193,]    0.4021   2.4870   17.5899
[194,]    0.4042   2.4742   17.7111
[195,]    0.4062   2.4615   17.4108
[196,]    0.4083   2.4490   17.2912
[197,]    0.4104   2.4365   16.7260
[198,]    0.4125   2.4242   16.3159
[199,]    0.4146   2.4121   15.8501
[200,]    0.4167   2.4000   15.4962
[201,]    0.4188   2.3881   14.6751
[202,]    0.4208   2.3762   13.9146
[203,]    0.4229   2.3645   13.1885
[204,]    0.4250   2.3529   13.2600
[205,]    0.4271   2.3415   12.7958
[206,]    0.4292   2.3301   12.5930
[207,]    0.4312   2.3188   12.8403
[208,]    0.4333   2.3077   13.1415
[209,]    0.4354   2.2967   13.2417
[210,]    0.4375   2.2857   14.0035
[211,]    0.4396   2.2749   14.8994
[212,]    0.4417   2.2642   15.9015
[213,]    0.4438   2.2535   16.0433
[214,]    0.4458   2.2430   16.5519
[215,]    0.4479   2.2326   17.0329
[216,]    0.4500   2.2222   16.0607
[217,]    0.4521   2.2120   15.3824
[218,]    0.4542   2.2018   14.8945
[219,]    0.4562   2.1918   13.8862
[220,]    0.4583   2.1818   13.0781
[221,]    0.4604   2.1719   12.3491
[222,]    0.4625   2.1622   12.0749
[223,]    0.4646   2.1525   11.9012
[224,]    0.4667   2.1429   11.7565
[225,]    0.4688   2.1333   12.5923
[226,]    0.4708   2.1239   13.1748
[227,]    0.4729   2.1145   13.5462
[228,]    0.4750   2.1053   14.0356
[229,]    0.4771   2.0961   14.2253
[230,]    0.4792   2.0870   14.0273
[231,]    0.4812   2.0779   13.6965
[232,]    0.4833   2.0690   13.1697
[233,]    0.4854   2.0601   12.5730
[234,]    0.4875   2.0513   11.8617
[235,]    0.4896   2.0426   11.0409
[236,]    0.4917   2.0339   10.3763
[237,]    0.4937   2.0253    9.7607
[238,]    0.4958   2.0168    9.3193
[239,]    0.4979   2.0084    9.1320
[240,]    0.5000   2.0000    9.0994

Cluster analysis


mydata <- read.csv("/Users/ansonleung/Desktop/Applied Financial Engineering/Financial Times Series Analysis/Lecture 3/utilities.csv")
head(mydata)

str(mydata)
'data.frame':   22 obs. of  9 variables:
 $ Company     : chr  "Arizona " "Boston " "Central " "Commonwealth" ...
 $ Fixed_charge: num  1.06 0.89 1.43 1.02 1.49 1.32 1.22 1.1 1.34 1.12 ...
 $ RoR         : num  9.2 10.3 15.4 11.2 8.8 13.5 12.2 9.2 13 12.4 ...
 $ Cost        : int  151 202 113 168 192 111 175 245 168 197 ...
 $ Load        : num  54.4 57.9 53 56 51.2 60 67.6 57 60.4 53 ...
 $ D.Demand    : num  1.6 2.2 3.4 0.3 1 -2.2 2.2 3.3 7.2 2.7 ...
 $ Sales       : int  9077 5088 9212 6423 3300 11127 7642 13082 8406 6455 ...
 $ Nuclear     : num  0 25.3 0 34.3 15.6 22.5 0 0 0 39.2 ...
 $ Fuel_Cost   : num  0.628 1.555 1.058 0.7 2.044 ...
pairs(mydata[,-c(1,1)])


plot(mydata$Fuel_Cost~ mydata$Sales, data = mydata)

with(mydata,
     text(mydata$Fuel_Cost ~ mydata$Sales, labels=mydata$Company,pos=4, cex = 0.4)
     )


z = mydata[,-c(1,1)]

means = apply(z,2,mean)
sds = apply(z,2,sd)
nor = scale(z,center=means,scale=sds)

distance = dist(nor)

print(distance, digits = 3)
      1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20   21
2  3.10                                                                                                    
3  3.68 4.92                                                                                               
4  2.46 2.16 4.11                                                                                          
5  4.12 3.85 4.47 4.13                                                                                     
6  3.61 4.22 2.99 3.20 4.60                                                                                
7  3.90 3.45 4.22 3.97 4.60 3.35                                                                           
8  2.74 3.89 4.99 3.69 5.16 4.91 4.36                                                                      
9  3.25 3.96 2.75 3.75 4.49 3.73 2.80 3.59                                                                 
10 3.10 2.71 3.93 1.49 4.05 3.83 4.51 3.67 3.57                                                            
11 3.49 4.79 5.90 4.86 6.46 6.00 6.00 3.46 5.18 5.08                                                       
12 3.22 2.43 4.03 3.50 3.60 3.74 1.66 4.06 2.74 3.94 5.21                                                  
13 3.96 3.43 4.39 2.58 4.76 4.55 5.01 4.14 3.66 1.41 5.31 4.50                                             
14 2.11 4.32 2.74 3.23 4.82 3.47 4.91 4.34 3.82 3.61 4.32 4.34 4.39                                        
15 2.59 2.50 5.16 3.19 4.26 4.07 2.93 3.85 4.11 4.26 4.74 2.33 5.10 4.24                                   
16 4.03 4.84 5.26 4.97 5.82 5.84 5.04 2.20 3.63 4.53 3.43 4.62 4.41 5.17 5.18                              
17 4.40 3.62 6.36 4.89 5.63 6.10 4.58 5.43 4.90 5.48 4.75 3.50 5.61 5.56 3.40 5.56                         
18 1.88 2.90 2.72 2.65 4.34 2.85 2.95 3.24 2.43 3.07 3.95 2.45 3.78 2.30 3.00 3.97 4.43                    
19 2.41 4.63 3.18 3.46 5.13 2.58 4.52 4.11 4.11 4.13 4.52 4.41 5.01 1.88 4.03 5.23 6.09 2.47               
20 3.17 3.00 3.73 1.82 4.39 2.91 3.54 4.09 2.95 2.05 5.35 3.43 2.23 3.74 3.78 4.82 4.87 2.92 3.90          
21 3.45 2.32 5.09 3.88 3.64 4.63 2.68 3.98 3.74 4.36 4.88 1.38 4.94 4.93 2.10 4.57 3.10 3.19 4.97 4.15     
22 2.51 2.42 4.11 2.58 3.77 4.03 4.00 3.24 3.21 2.56 3.44 3.00 2.74 3.51 3.35 3.46 3.63 2.55 3.97 2.62 3.01
mydata.hclust = hclust(distance)
plot(mydata.hclust)


plot(mydata.hclust,labels=mydata$Company,main='Default from hclust')


plot(mydata.hclust,hang=-1)


member = cutree(mydata.hclust,3)

table(member)
member
 1  2  3 
14  5  3 
aggregate(nor,list(member),mean)


aggregate(mydata[,-c(1,1)],list(member),mean)


wss <- (nrow(nor)-1)*sum(apply(nor,2,var))
for (i in 2:20) wss[i] <- sum(kmeans(nor, centers=i)$withinss)
plot(1:20, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")


kc <- kmeans(z,3)
plot(Sales~D.Demand, mydata, col = kc$cluster)

fractional R: Fractional Differences:

varve = scan("/Users/ansonleung/Desktop/Applied Financial Engineering/Financial Times Series Analysis/Lecture 3/varve.dat")
Read 634 items
varve = ts(varve)

library(arfima)

y = log(varve) - mean(log(varve))

acf(y)


varvefd = arfima(y)
Note: only one starting point.  Only one mode can be found -- this is now the default behavior.
Beginning the fits with 1 starting values.
d = summary(varvefd)$coef[[1]][1]

resids = resid(varvefd)[[1]] 

plot.ts(resids)

acf(resids)

Looping through p and q


data(Mishkin,package="Ecdat")

x <- diff(as.vector(Mishkin[,1]))

plot(x, type = 'l')


result <- matrix(0, nrow=9, ncol=4)
idx <- 1

for (i in 0:2)
{
  for (j in 0:2){
    fit = arima(x,order=c(i,0,j))
    result[idx, 1] = i
    result[idx, 2] = j
    result[idx, 3] = fit$aic
    result[idx, 4] = result[idx,3] + (log(length(x))-2)*i
    idx = idx + 1
    }
}
result <- data.frame(result)
names(result) <- c('p', 'q', 'AIC', 'BIC')
result

tar model


flu = scan("/Users/ansonleung/Desktop/Applied Financial Engineering/Financial Times Series Analysis/Lecture 3/flu.dat")
Read 132 items
flu = ts(flu)

plot(flu,type="b")


y = diff(flu,1)
plot(y,type="b")


model = ts.intersect(y, lag1y=lag(y,-1), lag2y=lag(y, -2), lag3y=lag(y,-3), lag4y=lag(y, -4))

x = model[,1]
P = model[,2:5]
c = 0.05 ## Threshold value

less = (P[,1]<c)
x1 = x[less]
P1 = P[less,]
out1 = lm(x1~P1[,1]+P1[,2]+P1[,3]+P1[,4])
summary(out1)

Call:
lm(formula = x1 ~ P1[, 1] + P1[, 2] + P1[, 3] + P1[, 4])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.13312 -0.02049  0.00218  0.01667  0.26666 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.004471   0.004894   0.914 0.363032    
P1[, 1]      0.506650   0.078319   6.469  3.2e-09 ***
P1[, 2]     -0.200086   0.056573  -3.537 0.000604 ***
P1[, 3]      0.121047   0.054463   2.223 0.028389 *  
P1[, 4]     -0.110938   0.045979  -2.413 0.017564 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04578 on 105 degrees of freedom
Multiple R-squared:  0.3763,    Adjusted R-squared:  0.3526 
F-statistic: 15.84 on 4 and 105 DF,  p-value: 3.568e-10
greater = (P[,1]>=c)
x2 = x[greater]
P2 = P[greater,]
out2 = lm(x2~P2[,1]+P2[,2]+P2[,3]+P2[,4])
summary(out2)

Call:
lm(formula = x2 ~ P2[, 1] + P2[, 2] + P2[, 3] + P2[, 4])

Residuals:
      Min        1Q    Median        3Q       Max 
-0.089975 -0.036825 -0.006328  0.040765  0.129509 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.40794    0.04675   8.726 1.53e-06 ***
P2[, 1]     -0.74833    0.16644  -4.496 0.000732 ***
P2[, 2]     -1.03231    0.21137  -4.884 0.000376 ***
P2[, 3]     -2.04504    1.05000  -1.948 0.075235 .  
P2[, 4]     -6.71178    1.24538  -5.389 0.000163 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0721 on 12 degrees of freedom
Multiple R-squared:  0.9207,    Adjusted R-squared:  0.8943 
F-statistic: 34.85 on 4 and 12 DF,  p-value: 1.618e-06
res1 = residuals(out1)
res2 = residuals(out2)
less[less==1] = res1
greater[greater==1] = res2
resid = less + greater
acf(resid)


less = (P[,1]<c)
greater = (P[,1]>=c)
fit1 = predict(out1)
fit2 = predict(out2)
less[less==1]= fit1
greater[greater==1] = fit2
fit = less + greater
plot(y, type="o")
lines(fit, col = "red", lty="dashed")

NA
NA


library(tsDyn)
dflu = diff(flu,1)
flu.tar4.05 = setar(dflu, m=4, thDelay=0, th=.05)

With the threshold you gave (0.05) there is a regime with less than trim=15% observations (86.61%, 13.39%, )

 1 T: Trim not respected:  0.8661417 0.1338583 from th: 0.05
Possible unit root in the high  regime. Roots are: 0.6182 0.6244 0.6244 0.6182
summary(flu.tar4.05)

Non linear autoregressive model

SETAR model ( 2 regimes)
Coefficients:
Low regime:
     const.L       phiL.1       phiL.2       phiL.3       phiL.4 
 0.004471044  0.506649694 -0.200086031  0.121047354 -0.110938271 

High regime:
   const.H     phiH.1     phiH.2     phiH.3     phiH.4 
 0.4079353 -0.7483325 -1.0323129 -2.0450407 -6.7117769 

Threshold:
-Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
-Value: 0.05 (fixed)
Proportion of points in low regime: 86.61%   High regime: 13.39% 

Residuals:
       Min         1Q     Median         3Q        Max 
-0.1331155 -0.0215013  0.0021557  0.0197805  0.2666619 

Fit:
residuals variance = 0.002156,  AIC = -784, MAPE = 404.3%

Coefficient(s):

          Estimate  Std. Error  t value  Pr(>|t|)    
const.L  0.0044710   0.0051648   0.8657  0.388381    
phiL.1   0.5066497   0.0826526   6.1299 1.137e-08 ***
phiL.2  -0.2000860   0.0597035  -3.3513  0.001073 ** 
phiL.3   0.1210474   0.0574764   2.1060  0.037268 *  
phiL.4  -0.1109383   0.0485235  -2.2863  0.023977 *  
const.H  0.4079353   0.0313281  13.0214 < 2.2e-16 ***
phiH.1  -0.7483325   0.1115341  -6.7094 6.658e-10 ***
phiH.2  -1.0323129   0.1416409  -7.2882 3.525e-11 ***
phiH.3  -2.0450407   0.7036312  -2.9064  0.004349 ** 
phiH.4  -6.7117769   0.8345586  -8.0423 6.771e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)+ (0) X(t-3)

Value: 0.05 (fixed)
plot(flu.tar4.05)

NA

NA

NA

NA

NA

NA

NA

NA

flu.tar4 = setar(dflu, m=4, thDelay=0)
Possible unit root in the high  regime. Roots are: 0.8074 0.8074 1.3596 1.3596
summary(flu.tar4)

Non linear autoregressive model

SETAR model ( 2 regimes)
Coefficients:
Low regime:
      const.L        phiL.1        phiL.2        phiL.3        phiL.4 
 0.0006269563  0.4608089284 -0.2243720404  0.1100931813 -0.1307031988 

High regime:
   const.H     phiH.1     phiH.2     phiH.3     phiH.4 
 0.2035231 -0.4071318 -1.4686776  0.3768388 -0.8298225 

Threshold:
-Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
-Value: 0.03646
Proportion of points in low regime: 84.25%   High regime: 15.75% 

Residuals:
       Min         1Q     Median         3Q        Max 
-0.2632248 -0.0195074  0.0044446  0.0210898  0.2854727 

Fit:
residuals variance = 0.003739,  AIC = -710, MAPE = 541.4%

Coefficient(s):

           Estimate  Std. Error  t value  Pr(>|t|)    
const.L  0.00062696  0.00699184   0.0897  0.928698    
phiL.1   0.46080893  0.11083474   4.1576 6.035e-05 ***
phiL.2  -0.22437204  0.07882706  -2.8464  0.005196 ** 
phiL.3   0.11009318  0.07608021   1.4471  0.150464    
phiL.4  -0.13070320  0.06406509  -2.0402  0.043511 *  
const.H  0.20352307  0.02632410   7.7314 3.505e-12 ***
phiH.1  -0.40713175  0.13413815  -3.0352  0.002944 ** 
phiH.2  -1.46867761  0.17439989  -8.4213 8.911e-14 ***
phiH.3   0.37683878  0.75269827   0.5007  0.617527    
phiH.4  -0.82982249  0.84051869  -0.9873  0.325478    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold
Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)+ (0) X(t-3)

Value: 0.03646
plot(flu.tar4)

NA

NA

NA

NA

NA

NA

NA

NA

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhpcyBpcyBhbiBbUiBNYXJrZG93bl0oaHR0cDovL3JtYXJrZG93bi5yc3R1ZGlvLmNvbSkgTm90ZWJvb2suIFdoZW4geW91IGV4ZWN1dGUgY29kZSB3aXRoaW4gdGhlIG5vdGVib29rLCB0aGUgcmVzdWx0cyBhcHBlYXIgYmVuZWF0aCB0aGUgY29kZS4gCgpUcnkgZXhlY3V0aW5nIHRoaXMgY2h1bmsgYnkgY2xpY2tpbmcgdGhlICpSdW4qIGJ1dHRvbiB3aXRoaW4gdGhlIGNodW5rIG9yIGJ5IHBsYWNpbmcgeW91ciBjdXJzb3IgaW5zaWRlIGl0IGFuZCBwcmVzc2luZyAqQ21kK1NoaWZ0K0VudGVyKi4gCgoKQWRkIGEgbmV3IGNodW5rIGJ5IGNsaWNraW5nIHRoZSAqSW5zZXJ0IENodW5rKiBidXR0b24gb24gdGhlIHRvb2xiYXIgb3IgYnkgcHJlc3NpbmcgKkNtZCtPcHRpb24rSSouCgpXaGVuIHlvdSBzYXZlIHRoZSBub3RlYm9vaywgYW4gSFRNTCBmaWxlIGNvbnRhaW5pbmcgdGhlIGNvZGUgYW5kIG91dHB1dCB3aWxsIGJlIHNhdmVkIGFsb25nc2lkZSBpdCAoY2xpY2sgdGhlICpQcmV2aWV3KiBidXR0b24gb3IgcHJlc3MgKkNtZCtTaGlmdCtLKiB0byBwcmV2aWV3IHRoZSBIVE1MIGZpbGUpLiAKClRoZSBwcmV2aWV3IHNob3dzIHlvdSBhIHJlbmRlcmVkIEhUTUwgY29weSBvZiB0aGUgY29udGVudHMgb2YgdGhlIGVkaXRvci4gQ29uc2VxdWVudGx5LCB1bmxpa2UgKktuaXQqLCAqUHJldmlldyogZG9lcyBub3QgcnVuIGFueSBSIGNvZGUgY2h1bmtzLiBJbnN0ZWFkLCB0aGUgb3V0cHV0IG9mIHRoZSBjaHVuayB3aGVuIGl0IHdhcyBsYXN0IHJ1biBpbiB0aGUgZWRpdG9yIGlzIGRpc3BsYXllZC4KCkstbWVhbgpgYGB7cn0KaGVhZChpcmlzKQpwbG90KGlyaXNbLDM6NF0pCgppcmlzQ2x1c3RlciA8LSBrbWVhbnMoaXJpc1ssIDM6NF0sIDMsIG5zdGFydCA9IDIwKQppcmlzQ2x1c3RlciRjbHVzdGVyIDwtIGFzLmZhY3RvcihpcmlzQ2x1c3RlciRjbHVzdGVyKQpsaWJyYXJ5KGdncGxvdDIpCgpnZ3Bsb3QoaXJpcywgYWVzKFBldGFsLkxlbmd0aCwgUGV0YWwuV2lkdGgsIGNvbCA9IGlyaXNDbHVzdGVyJGNsdXN0ZXIsIGZpbGw9U3BlY2llcykpICsgCiAgc3RhdF9lbGxpcHNlKGdlb20gPSAicG9seWdvbiIsIGNvbCA9ICJibGFjayIsIGFscGhhID0gMC41KSArCiAgZ2VvbV9wb2ludCgpCgprLm1heCA8LSAxNQp3c3MgPC0gc2FwcGx5KDE6ay5tYXgsIGZ1bmN0aW9uKGspe2ttZWFucyhpcmlzWywgMzo0XSwgaywgbnN0YXJ0PTUwLCBpdGVyLm1heCA9IDE1KSR0b3Qud2l0aGluc3N9KQoKcGxvdCgxOmsubWF4LCB3c3MsCiAgICAgdHlwZT0iYiIsIHBjaCA9IDE5LCBmcmFtZSA9IEZBTFNFLCAKICAgICB4bGFiPSJOdW1iZXIgb2YgY2x1c3RlcnMgSyIsCiAgICAgeWxhYj0iVG90YWwgd2l0aGluLWNsdXN0ZXJzIHN1bSBvZiBzcXVhcmVzIikKCmBgYAoKCgpQZXJpb2RvZ3JhbQpgYGB7cn0KY291bnRlciA8LSA1MDAKeCA8LSBzZXEoY291bnRlcikKCnkgPC0gMiAqIGNvcygyICogcGkgKiAoMS81MCkgKiB4ICsgMC42ICogcGkpCnNpZ25hbCA8LSB5ICsgcm5vcm0oY291bnRlcikKcGxvdCh4LCB5LCB0eXBlID0gImwiKQoKcGxvdCh4LCBzaWduYWwsIHR5cGUgPSAibCIpCmBgYApgYGB7cn0Kc3Vuc3BvdHMgPC0gc2NhbigiL1VzZXJzL2Fuc29ubGV1bmcvRGVza3RvcC9BcHBsaWVkIEZpbmFuY2lhbCBFbmdpbmVlcmluZy9GaW5hbmNpYWwgVGltZXMgU2VyaWVzIEFuYWx5c2lzL0xlY3R1cmUgMy9zdW5zcG90cy5kYXQiKQpwbG90KHN1bnNwb3RzLCB0eXBlPSJiIikKeCA9IGRpZmYoc3Vuc3BvdHMpCgpJID0gYWJzKGZmdCh4KS9zcXJ0KDQ1OCkpXjIKUCA9ICg0LzQ1OCkqSVsxOjIzMF0KCmZyZXEgPSAoMDoyMjkpLzQ1OAoKcGxvdChmcmVxLCBQLCB0eXBlPSJsIikKYGBgCgoKCmBgYHtyfQoKbGlicmFyeShhc3RzYSkKCnggPC0gc2NhbigiL1VzZXJzL2Fuc29ubGV1bmcvRGVza3RvcC9BcHBsaWVkIEZpbmFuY2lhbCBFbmdpbmVlcmluZy9GaW5hbmNpYWwgVGltZXMgU2VyaWVzIEFuYWx5c2lzL0xlY3R1cmUgMy9yZWNydWl0LmRhdCIpCm12c3BlYyh4LCBsb2c9Im5vIikKCmsgPSBrZXJuZWwoImRhbmllbGwiLCA0KQptdnNwZWMoeCwgaywgbG9nPSJubyIpCgprID0ga2VybmVsKCJkYW5pZWxsIiwgYyg0LDQpKSAKbXZzcGVjKHgsIGssIGxvZz0ibm8iKQoKCnNwZWN2YWx1ZXMgPSBtdnNwZWMoeCwgaywgbG9nPSJubyIpCnNwZWN2YWx1ZXMkZGV0YWlscwoKYGBgCgpDbHVzdGVyIGFuYWx5c2lzCmBgYHtyfQoKbXlkYXRhIDwtIHJlYWQuY3N2KCIvVXNlcnMvYW5zb25sZXVuZy9EZXNrdG9wL0FwcGxpZWQgRmluYW5jaWFsIEVuZ2luZWVyaW5nL0ZpbmFuY2lhbCBUaW1lcyBTZXJpZXMgQW5hbHlzaXMvTGVjdHVyZSAzL3V0aWxpdGllcy5jc3YiKQpoZWFkKG15ZGF0YSkKCnN0cihteWRhdGEpCnBhaXJzKG15ZGF0YVssLWMoMSwxKV0pCgpwbG90KG15ZGF0YSRGdWVsX0Nvc3R+IG15ZGF0YSRTYWxlcywgZGF0YSA9IG15ZGF0YSkKCndpdGgobXlkYXRhLAogICAgIHRleHQobXlkYXRhJEZ1ZWxfQ29zdCB+IG15ZGF0YSRTYWxlcywgbGFiZWxzPW15ZGF0YSRDb21wYW55LHBvcz00LCBjZXggPSAwLjQpCiAgICAgKQoKeiA9IG15ZGF0YVssLWMoMSwxKV0KCm1lYW5zID0gYXBwbHkoeiwyLG1lYW4pCnNkcyA9IGFwcGx5KHosMixzZCkKbm9yID0gc2NhbGUoeixjZW50ZXI9bWVhbnMsc2NhbGU9c2RzKQoKZGlzdGFuY2UgPSBkaXN0KG5vcikKCnByaW50KGRpc3RhbmNlLCBkaWdpdHMgPSAzKQoKbXlkYXRhLmhjbHVzdCA9IGhjbHVzdChkaXN0YW5jZSkKcGxvdChteWRhdGEuaGNsdXN0KQoKcGxvdChteWRhdGEuaGNsdXN0LGxhYmVscz1teWRhdGEkQ29tcGFueSxtYWluPSdEZWZhdWx0IGZyb20gaGNsdXN0JykKCnBsb3QobXlkYXRhLmhjbHVzdCxoYW5nPS0xKQoKbWVtYmVyID0gY3V0cmVlKG15ZGF0YS5oY2x1c3QsMykKCnRhYmxlKG1lbWJlcikKCmFnZ3JlZ2F0ZShub3IsbGlzdChtZW1iZXIpLG1lYW4pCgoKYWdncmVnYXRlKG15ZGF0YVssLWMoMSwxKV0sbGlzdChtZW1iZXIpLG1lYW4pCgoKd3NzIDwtIChucm93KG5vciktMSkqc3VtKGFwcGx5KG5vciwyLHZhcikpCmZvciAoaSBpbiAyOjIwKSB3c3NbaV0gPC0gc3VtKGttZWFucyhub3IsIGNlbnRlcnM9aSkkd2l0aGluc3MpCnBsb3QoMToyMCwgd3NzLCB0eXBlPSJiIiwgeGxhYj0iTnVtYmVyIG9mIENsdXN0ZXJzIiwgeWxhYj0iV2l0aGluIGdyb3VwcyBzdW0gb2Ygc3F1YXJlcyIpCgprYyA8LSBrbWVhbnMoeiwzKQpwbG90KFNhbGVzfkQuRGVtYW5kLCBteWRhdGEsIGNvbCA9IGtjJGNsdXN0ZXIpCgpgYGAKCmZyYWN0aW9uYWwgUjogRnJhY3Rpb25hbCBEaWZmZXJlbmNlczoKYGBge3J9CnZhcnZlID0gc2NhbigiL1VzZXJzL2Fuc29ubGV1bmcvRGVza3RvcC9BcHBsaWVkIEZpbmFuY2lhbCBFbmdpbmVlcmluZy9GaW5hbmNpYWwgVGltZXMgU2VyaWVzIEFuYWx5c2lzL0xlY3R1cmUgMy92YXJ2ZS5kYXQiKQp2YXJ2ZSA9IHRzKHZhcnZlKQoKbGlicmFyeShhcmZpbWEpCgp5ID0gbG9nKHZhcnZlKSAtIG1lYW4obG9nKHZhcnZlKSkKCmFjZih5KQoKdmFydmVmZCA9IGFyZmltYSh5KQoKZCA9IHN1bW1hcnkodmFydmVmZCkkY29lZltbMV1dWzFdCgpyZXNpZHMgPSByZXNpZCh2YXJ2ZWZkKVtbMV1dIAoKcGxvdC50cyhyZXNpZHMpCmFjZihyZXNpZHMpCmBgYAoKTG9vcGluZyB0aHJvdWdoIHAgYW5kIHEKYGBge3J9CgpkYXRhKE1pc2hraW4scGFja2FnZT0iRWNkYXQiKQoKeCA8LSBkaWZmKGFzLnZlY3RvcihNaXNoa2luWywxXSkpCgpwbG90KHgsIHR5cGUgPSAnbCcpCgpyZXN1bHQgPC0gbWF0cml4KDAsIG5yb3c9OSwgbmNvbD00KQppZHggPC0gMQoKZm9yIChpIGluIDA6MikKewogIGZvciAoaiBpbiAwOjIpewogICAgZml0ID0gYXJpbWEoeCxvcmRlcj1jKGksMCxqKSkKICAgIHJlc3VsdFtpZHgsIDFdID0gaQogICAgcmVzdWx0W2lkeCwgMl0gPSBqCiAgICByZXN1bHRbaWR4LCAzXSA9IGZpdCRhaWMKICAgIHJlc3VsdFtpZHgsIDRdID0gcmVzdWx0W2lkeCwzXSArIChsb2cobGVuZ3RoKHgpKS0yKSppCiAgICBpZHggPSBpZHggKyAxCiAgICB9Cn0KcmVzdWx0IDwtIGRhdGEuZnJhbWUocmVzdWx0KQpuYW1lcyhyZXN1bHQpIDwtIGMoJ3AnLCAncScsICdBSUMnLCAnQklDJykKcmVzdWx0CmBgYAoKCnRhciBtb2RlbApgYGB7cn0KCmZsdSA9IHNjYW4oIi9Vc2Vycy9hbnNvbmxldW5nL0Rlc2t0b3AvQXBwbGllZCBGaW5hbmNpYWwgRW5naW5lZXJpbmcvRmluYW5jaWFsIFRpbWVzIFNlcmllcyBBbmFseXNpcy9MZWN0dXJlIDMvZmx1LmRhdCIpCmZsdSA9IHRzKGZsdSkKCnBsb3QoZmx1LHR5cGU9ImIiKQoKeSA9IGRpZmYoZmx1LDEpCnBsb3QoeSx0eXBlPSJiIikKCm1vZGVsID0gdHMuaW50ZXJzZWN0KHksIGxhZzF5PWxhZyh5LC0xKSwgbGFnMnk9bGFnKHksIC0yKSwgbGFnM3k9bGFnKHksLTMpLCBsYWc0eT1sYWcoeSwgLTQpKQoKeCA9IG1vZGVsWywxXQpQID0gbW9kZWxbLDI6NV0KYyA9IDAuMDUgIyMgVGhyZXNob2xkIHZhbHVlCgpsZXNzID0gKFBbLDFdPGMpCngxID0geFtsZXNzXQpQMSA9IFBbbGVzcyxdCm91dDEgPSBsbSh4MX5QMVssMV0rUDFbLDJdK1AxWywzXStQMVssNF0pCnN1bW1hcnkob3V0MSkKCmdyZWF0ZXIgPSAoUFssMV0+PWMpCngyID0geFtncmVhdGVyXQpQMiA9IFBbZ3JlYXRlcixdCm91dDIgPSBsbSh4Mn5QMlssMV0rUDJbLDJdK1AyWywzXStQMlssNF0pCnN1bW1hcnkob3V0MikKCnJlczEgPSByZXNpZHVhbHMob3V0MSkKcmVzMiA9IHJlc2lkdWFscyhvdXQyKQpsZXNzW2xlc3M9PTFdID0gcmVzMQpncmVhdGVyW2dyZWF0ZXI9PTFdID0gcmVzMgpyZXNpZCA9IGxlc3MgKyBncmVhdGVyCmFjZihyZXNpZCkKCmxlc3MgPSAoUFssMV08YykKZ3JlYXRlciA9IChQWywxXT49YykKZml0MSA9IHByZWRpY3Qob3V0MSkKZml0MiA9IHByZWRpY3Qob3V0MikKbGVzc1tsZXNzPT0xXT0gZml0MQpncmVhdGVyW2dyZWF0ZXI9PTFdID0gZml0MgpmaXQgPSBsZXNzICsgZ3JlYXRlcgpwbG90KHksIHR5cGU9Im8iKQpsaW5lcyhmaXQsIGNvbCA9ICJyZWQiLCBsdHk9ImRhc2hlZCIpCgoKYGBgCgoKYGBge3J9CgoKbGlicmFyeSh0c0R5bikKZGZsdSA9IGRpZmYoZmx1LDEpCmZsdS50YXI0LjA1ID0gc2V0YXIoZGZsdSwgbT00LCB0aERlbGF5PTAsIHRoPS4wNSkKc3VtbWFyeShmbHUudGFyNC4wNSkKCnBsb3QoZmx1LnRhcjQuMDUpCgpgYGAKCgpgYGB7cn0KZmx1LnRhcjQgPSBzZXRhcihkZmx1LCBtPTQsIHRoRGVsYXk9MCkKc3VtbWFyeShmbHUudGFyNCkKcGxvdChmbHUudGFyNCkKYGBgCgpgYGB7cn0KCgoKYGBgCgoKYGBge3J9CgoKCgoKYGBgCgoKCgo=